home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / LLSQ.FOR < prev    next >
Text File  |  1985-11-29  |  6KB  |  238 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE LLSQ
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO MINIMIZE
  8. C           THE EUCLIDEAN NORM OF B-A*X, WHERE A IS A M BY N MATRIX
  9. C           WITH M NOT LESS THAN N. IN THE SPECIAL CASE M=N SYSTEMS OF
  10. C           LINEAR EQUATIONS MAY BE SOLVED.
  11. C
  12. C        USAGE
  13. C           CALL LLSQ (A,B,M,N,L,X,IPIV,EPS,IER,AUX)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A      - M BY N COEFFICIENT MATRIX (DESTROYED).
  17. C           B      - M BY L RIGHT HAND SIDE MATRIX (DESTROYED).
  18. C           M      - ROW NUMBER OF MATRICES A AND B.
  19. C           N      - COLUMN NUMBER OF MATRIX A, ROW NUMBER OF MATRIX X.
  20. C           L      - COLUMN NUMBER OF MATRICES B AND X.
  21. C           X      - N BY L SOLUTION MATRIX.
  22. C           IPIV   - INTEGER OUTPUT VECTOR OF DIMENSION N WHICH
  23. C                    CONTAINS INFORMATIONS ON COLUMN INTERCHANGES
  24. C                    IN MATRIX A. (SEE REMARK NO.3).
  25. C           EPS    - INPUT PARAMETER WHICH SPECIFIES A RELATIVE
  26. C                    TOLERANCE FOR DETERMINATION OF RANK OF MATRIX A.
  27. C           IER    - A RESULTING ERROR PARAMETER.
  28. C           AUX    - AUXILIARY STORAGE ARRAY OF DIMENSION MAX(2*N,L).
  29. C                    ON RETURN FIRST L LOCATIONS OF AUX CONTAIN THE
  30. C                    RESULTING LEAST SQUARES.
  31. C
  32. C        REMARKS
  33. C           (1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN CASE
  34. C               M LESS THAN N.
  35. C           (2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN CASE
  36. C               OF A ZERO-MATRIX A.
  37. C           (3) IF RANK K OF MATRIX A IS FOUND TO BE LESS THAN N BUT
  38. C               GREATER THAN 0, THE PROCEDURE RETURNS WITH ERROR CODE
  39. C               IER=K INTO CALLING PROGRAM. THE LAST N-K ELEMENTS OF
  40. C               VECTOR IPIV DENOTE THE USELESS COLUMNS IN MATRIX A.
  41. C               THE REMAINING USEFUL COLUMNS FORM A BASE OF MATRIX A.
  42. C           (4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR PARAMETER IER
  43. C               IS SET TO 0.
  44. C
  45. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  46. C           NONE
  47. C
  48. C        METHOD
  49. C           HOUSEHOLDER TRANSFORMATIONS ARE USED TO TRANSFORM MATRIX A
  50. C           TO UPPER TRIANGULAR FORM. AFTER HAVING APPLIED THE SAME
  51. C           TRANSFORMATION TO THE RIGHT HAND SIDE MATRIX B, AN
  52. C           APPROXIMATE SOLUTION OF THE PROBLEM IS COMPUTED BY
  53. C           BACK SUBSTITUTION. FOR REFERENCE, SEE
  54. C           G. GOLUB, NUMERICAL METHODS FOR SOLVING LINEAR LEAST
  55. C           SQUARES PROBLEMS, NUMERISCHE MATHEMATIK, VOL.7,
  56. C           ISS.3 (1965), PP.206-216.
  57. C
  58. C     ..................................................................
  59. C
  60.       SUBROUTINE LLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)
  61. C
  62.       DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)
  63. C
  64. C     ERROR TEST
  65.       IF(M-N)30,1,1
  66. C
  67. C     GENERATION OF INITIAL VECTOR S(K) (K=1,2,...,N) IN STORAGE
  68. C     LOCATIONS AUX(K) (K=1,2,...,N)
  69.     1 PIV=0.
  70.       IEND=0
  71.       DO 4 K=1,N
  72.       IPIV(K)=K
  73.       H=0.
  74.       IST=IEND+1
  75.       IEND=IEND+M
  76.       DO 2 I=IST,IEND
  77.     2 H=H+A(I)*A(I)
  78.       AUX(K)=H
  79.       IF(H-PIV)4,4,3
  80.     3 PIV=H
  81.       KPIV=K
  82.     4 CONTINUE
  83. C
  84. C     ERROR TEST
  85.       IF(PIV)31,31,5
  86. C
  87. C     DEFINE TOLERANCE FOR CHECKING RANK OF A
  88.     5 SIG=SQRT(PIV)
  89.       TOL=SIG*ABS(EPS)
  90. C
  91. C
  92. C     DECOMPOSITION LOOP
  93.       LM=L*M
  94.       IST=-M
  95.       DO 21 K=1,N
  96.       IST=IST+M+1
  97.       IEND=IST+M-K
  98.       I=KPIV-K
  99.       IF(I)8,8,6
  100. C
  101. C     INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K
  102.     6 H=AUX(K)
  103.       AUX(K)=AUX(KPIV)
  104.       AUX(KPIV)=H
  105.       ID=I*M
  106.       DO 7 I=IST,IEND
  107.       J=I+ID
  108.       H=A(I)
  109.       A(I)=A(J)
  110.     7 A(J)=H
  111. C
  112. C     COMPUTATION OF PARAMETER SIG
  113.     8 IF(K-1)11,11,9
  114.     9 SIG=0.
  115.       DO 10 I=IST,IEND
  116.    10 SIG=SIG+A(I)*A(I)
  117.       SIG=SQRT(SIG)
  118. C
  119. C     TEST ON SINGULARITY
  120.       IF(SIG-TOL)32,32,11
  121. C
  122. C     GENERATE CORRECT SIGN OF PARAMETER SIG
  123.    11 H=A(IST)
  124.       IF(H)12,13,13
  125.    12 SIG=-SIG
  126. C
  127. C     SAVE INTERCHANGE INFORMATION
  128.    13 IPIV(KPIV)=IPIV(K)
  129.       IPIV(K)=KPIV
  130. C
  131. C     GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF
  132. C     PARAMETER BETA
  133.       BETA=H+SIG
  134.       A(IST)=BETA
  135.       BETA=1./(SIG*BETA)
  136.       J=N+K
  137.       AUX(J)=-SIG
  138.       IF(K-N)14,19,19
  139. C
  140. C     TRANSFORMATION OF MATRIX A
  141.    14 PIV=0.
  142.       ID=0
  143.       JST=K+1
  144.       KPIV=JST
  145.       DO 18 J=JST,N
  146.       ID=ID+M
  147.       H=0.
  148.       DO 15 I=IST,IEND
  149.       II=I+ID
  150.    15 H=H+A(I)*A(II)
  151.       H=BETA*H
  152.       DO 16 I=IST,IEND
  153.       II=I+ID
  154.    16 A(II)=A(II)-A(I)*H
  155. C
  156. C     UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J)
  157.       II=IST+ID
  158.       H=AUX(J)-A(II)*A(II)
  159.       AUX(J)=H
  160.       IF(H-PIV)18,18,17
  161.    17 PIV=H
  162.       KPIV=J
  163.    18 CONTINUE
  164. C
  165. C     TRANSFORMATION OF RIGHT HAND SIDE MATRIX B
  166.    19 DO 21 J=K,LM,M
  167.       H=0.
  168.       IEND=J+M-K
  169.       II=IST
  170.       DO 20 I=J,IEND
  171.       H=H+A(II)*B(I)
  172.    20 II=II+1
  173.       H=BETA*H
  174.       II=IST
  175.       DO 21 I=J,IEND
  176.       B(I)=B(I)-A(II)*H
  177.    21 II=II+1
  178. C     END OF DECOMPOSITION LOOP
  179. C
  180. C
  181. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  182.       IER=0
  183.       I=N
  184.       LN=L*N
  185.       PIV=1./AUX(2*N)
  186.       DO 22 K=N,LN,N
  187.       X(K)=PIV*B(I)
  188.    22 I=I+M
  189.       IF(N-1)26,26,23
  190.    23 JST=(N-1)*M+N
  191.       DO 25 J=2,N
  192.       JST=JST-M-1
  193.       K=N+N+1-J
  194.       PIV=1./AUX(K)
  195.       KST=K-N
  196.       ID=IPIV(KST)-KST
  197.       IST=2-J
  198.       DO 25 K=1,L
  199.       H=B(KST)
  200.       IST=IST+N
  201.       IEND=IST+J-2
  202.       II=JST
  203.       DO 24 I=IST,IEND
  204.       II=II+M
  205.    24 H=H-A(II)*X(I)
  206.       I=IST-1
  207.       II=I+ID
  208.       X(I)=X(II)
  209.       X(II)=PIV*H
  210.    25 KST=KST+M
  211. C
  212. C
  213. C     COMPUTATION OF LEAST SQUARES
  214.    26 IST=N+1
  215.       IEND=0
  216.       DO 29 J=1,L
  217.       IEND=IEND+M
  218.       H=0.
  219.       IF(M-N)29,29,27
  220.    27 DO 28 I=IST,IEND
  221.    28 H=H+B(I)*B(I)
  222.       IST=IST+M
  223.    29 AUX(J)=H
  224.       RETURN
  225. C
  226. C     ERROR RETURN IN CASE M LESS THAN N
  227.    30 IER=-2
  228.       RETURN
  229. C
  230. C     ERROR RETURN IN CASE OF ZERO-MATRIX A
  231.    31 IER=-1
  232.       RETURN
  233. C
  234. C     ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N
  235.    32 IER=K-1
  236.       RETURN
  237.       END
  238.